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Abstract 

Theoretical frameworks to estimate the tolerance of metabolic networks to various 
failures are important to evaluate the robustness of biological complex systems in 
systems biology. In this paper, we focus on a measure for robustness in metabolic 
networks, namely, the impact degree, and propose an approximation method to pre- 
dict the probability distribution of impact degrees from metabolic network struc- 
tures using the theory of branching process. We demonstrate the relevance of this 
method by testing it on real-world metabolic networks. Although the approxima- 
tion method possesses a few limitations, it may be a powerful tool for evaluating 
metabolic robustness. 
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1 Introduction 



Robustness is a key feature in the analysis of complex systems, especially 
for complex biological systems. Many organisms have strong adaptability to 
environmental changes or failures in some of their components, and can live 
even if some of their genes are mutated. In particular, it is known that cancer 
cells are very robust pp. Therefore, understanding the origin of robustness of 
living cells has become an important research topic. 

In particular, extensive studies have focused on the analysis of structural ro- 
bustness of metabolic networks. Structural robustness refers to the tolerance 
of the system's behavior to changes in the structure of networks, and most 
existing studies focus on changes caused by knockout of gene(s) or enzyme(s). 
One of the reasons why extensive studies have been done on structural robust- 
ness of metabolic networks is that rather accurate data of metabolic networks 
are available via such databases as the Kyoto Encyclopedia of Genes and 
Genomes (KEGG) [2] and the Encyclopedia of Escherichia coli K-12 Genes 
and Metabolism (EcoCyc) [3J, and kinetic parameters, which are not neces- 
sarily available, are not required. 

In order to analyze the structural robustness of metabolic networks, the flux 
balance analysis (FBA) methods have been widely used. In many of these 
approaches, elementary flux modes (EFMs) play a key role, where an EFM 
is a minimal set of reactions that can operate at steady state [I]. Based on 
FBA and/or EFM, several studies have focused on finding a minimum reac- 
tion cut [5|6|7|8] , that is, a minimum set of reaction (or enzyme) removals 
which prevent the production of a specified set of compounds. Other FBA- 
based measures of robustness have also been proposed. Behre et al. proposed 
a measure based on the number of remaining EFMs after knockout versus the 
number of EFMs in the unperturbed situation [9J. Deutscher et al. proposed 
another measure using the Shaply value from game theory [ID] . 

Other approaches have been proposed based on Boolean models of metabolic 
networks in which reactions and compounds are modeled as AND and OR 
nodes, respectively. Handorf et al. analyzed robustness of metabolic networks 
by introducing the concept of scope [II]. Li et al., Sridhar et al., and Tamura 
et al. developed integer programming-based methods for finding a minimum 
reaction cut in Boolean models of metabolic networks [l2|13|14|fT5] . Smart et 
al. defined the topological flux balance (TFB) criterion based on a Boolean 
model of metabolic networks and analyzed the damage (number of reactions) 
caused by knockout of a single reaction under TFB [16]. Jiang et al. defined 
and analyzed the impact degree, which is the number of reactions inactivated 
by knockout of a specific reaction [IT] . Although there are some differences in 
the treatment of reversible reactions, the damage and the impact degree are 
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very similar concepts. Cong et al. extended the impact degree for knockout of 
multiple reactions [T8] . 

In this paper, we study the distribution of the impact degree caused by random 
knockout of a single reaction using the theory of branching process [T9|20] . In 
order to analyze earthquakes, Saichev et al. proposed a branching process 
with power-law distributions of offspring d: P{d) oc 1/gP +1 , where 7 is some 
constant, and approximately derived the distribution of the total number of 
offsprings [20]. We regard propagation of the impact of knockout of a reac- 
tion as a branching process, and apply their method to estimate the impact 
degree distribution, where the impact degree in our problem corresponds to 
the total number of offsprings in the branching process. In order to apply 
this method, we develop a simple method for estimating the offspring distri- 
bution in a metabolic network. Although Smart et al. have already applied 
percolation theory and branching process to analysis of the size distribution 
of rigid clusters, defined as clusters of contagion nodes that do not contain 
any branched metabolite nodes (see Fig. 3 A in [IS]), they did not explicitly 
estimate the damage distribution (i.e., the impact degree distribution). We 
finally show an estimation method for the damage distributions. The pro- 
posed method is applied to analysis of metabolic networks of four species: 
Escherichia coli, Bacillus subtilis, Saccharomyces cerevisiae and Homo sapi- 
ens. The results show good agreement of impact degree distributions between 
empirical results and theoretical estimates. 



2 Impact degree 

Jiang et al. proposed the impact degree as a measure of the importance of 
each reaction in a metabolic network [T7]. The impact degree is defined as 
the number of inactivated reactions caused by knockout of a single reaction. 
However, they did not consider the effect of cycles in metabolic networks. Since 
cycles play an important role in metabolic networks, Cong et al. extended the 
impact degree so that the effect of cycles is taken into account [18] by using 
a concept of the maximal valid assignment [15]. Here, we briefly review the 
definition of this extended impact degree [Tg] . 

As in other works, we regard each metabolic network as a bipartite directed 
graph. Let V c = {Ci, . . . , C m } and V r = {Ri, . . . , R n } be a set of compound 
nodes and a set of reaction nodes respectively, where V^flK- = {}. A metabolic 
network is defined as a directed graph G(V C U V r ,E) in which either (u G 
V c ) A (v G V r ) or (u G V r ) A (v G V c ) holds for each edge (u, v) G 
Each reaction and compound takes one of two states: or 1, where and 1 

1 A f\B means logical AND of A and B. 
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correspond to inactive and active reactions (compounds), respectively. Reverse 
reactions are treated as two irreversible reactions. 

In order to define the impact degree of a reaction, we proceed as follows. 
Suppose that reaction Ri is knocked out. Then, we start with the global state 
where all compounds are active (i.e., Ck = 1 for all e V c ) and all reactions 
but Ri are active (i.e., Rj = 1 for all Rj G V r \{Ri} and Ri = 0). Then, we 
alternatively update the states of reactions and compounds by the following 
rules. 

(1) For each reaction, there are three different compounds: consumed com- 
pounds (i.e., substrates), produced compounds (i.e., products), and di- 
rectly unrelated compounds. 

(2) A reaction is inactivated if any of its consumed or produced compound 
is inactivated. 

(3) For each compound, there are three different reactions: consuming reac- 
tions, producing reactions, and directly unrelated reaction. 

(4) A compound is inactivated if all its consuming reactions or all its pro- 
ducing reactions are inactivated. 

Since no activation is possible in this process, the procedure converges to a 
stable state in a finite number of iterations. The impact degree of reaction 
Ri is defined as the number of inactivated reactions in the stable state. This 
procedure simultaneously gives the definition of the impact degree and an 
algorithm to compute it. 




Fig. 1. Example of a metabolic network. Boxes and circles correspond to reactions 
and compounds, respectively. 

Let us illustrate the above process with the metabolic network shown in Fig. [Tl 
Suppose that reaction R 2 is knocked out. Then, the states of nodes change as 
shown in Table [U Since four reactions (including R 2 ) are inactivated in the 
stable state, the resulting impact degree is four. Next, suppose that reaction R e 
is knocked out. In this case, the states of nodes change as shown in Tableland 
the resulting impact degree is two. It is to be noted that C 4 is not inactivated 
because R 2 is still active in this case. 
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Table 1 

Impact degree calculation when R2 in Fig. Q] is knocked out. 
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Table 2 

Impact degree calculation when Rq in Fig. Q] is knocked out. 
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3 Branching process approximation 



We here explain the branching process approximation for estimating the im- 
pact degree distributions of metabolic networks. 

The branching process is a stochastic process in which each progenitor gener- 
ates offsprings according to a fixed probability distribution called the offspring 
distribution. We propose that the branching process approximation is useful 
for estimating the impact degree distribution because the propagation of an 
impact on a network is essentially similar to cascading failures, which is a 
sequence of failures caused by an accident. The branching process model of 
cascading failure is a standard Galton- Watson branching process [2T]. The 
approximation method using branching process has been already applied to 
loading-dependent cascading failure, and its relevance has been shown [22|23] . 
However, the branching process approximation needs the assumption of tree 
structure of networks; thus, this is a limitation of the branching process ap- 
proximation because metabolic networks generally have cycles. 
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3. 1 The number of offsprings in metabolic networks 



To estimate the impact degree distributions using the branching process ap- 
proximation, we need to define the notion of offsprings for each reaction node 
in metabolic networks, and the distribution of the number of offsprings. 

For that purpose, we consider the reaction network obtained as the unipartite 
projection of the metabolic network, where we draw an edge from reaction A 
to reaction B when at least one product of A is a substrate of B [25,26,2 7f28] . 
Fig. [2] shows the reaction network obtained from the metabolic network of Fig. 
[Qby this procedure. 

As an easy example, we consider 2 metabolic reactions, A and B. In this 
case, the edge is drawn from A to B (i.e., A— >B) if at least 1 product of 
reaction A corresponds to at least 1 substrate of reaction B (e.g., the case 
of a— ?• A— >b— >B— >c, where a, b, and c are chemical compounds). Through 
a similar procedure, we obtain a reaction network from a given metabolic 
network (see Fig. [2]). 




Fig. 2. The reaction network transformed from the metabolic network in Fig. [TJ 



We now observe that the number of reactions inactivated by the failure of a 
given reaction, which we want to model in the branching process, does not cor- 
respond to its outdegree (i.e., the number of out-going edges from a reaction 
node) in the reaction network because the spreading of an impact is depressed 
when there are alternative synthetic pathways. For example, assuming that 
the reaction Ri is inactive (or disrupted), the cascading of the impact does not 
occur because the reaction R4 remains active due to the chemical compound 
generated through the reaction R2. In contrast, the cascading of the impact 
continues when the reaction R2 becomes inactive for instance because the re- 
action Rq is dependent on this reaction only. In short, the impact spreads 
though reactions whose substrates are synthesized via unique metabolic reac- 
tions (i.e., reaction nodes with the indegree of 1) when assuming tree structures 
of networks. 



Based on this analysis, we define the number of offsprings for each reaction 
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node in metabolic networks as follows: 

Lout (jf fc in = ^ 

* = < (1) 
(otherwise) 



where k° ut and k\ n are the outdegree and indegree of reaction node % in the 
reaction network, respectively. 



3.2 Branching process models 



To analytically estimate the impact degree distributions, it is useful to assume 
that the number of offsprings for each node follows a probability distribution. 
We here consider two types of distributions, which are frequently observed in 
real worlds. 



3.2.1 Poisson model 

The simplest case of branching processes is a Poisson branching process (here- 
after called Poisson model) in which the number of offsprings d for each pro- 
genitor follows the Poisson distribution: /i d e~ M /<i!, where /x corresponds to the 
mean of this distribution. In this model, the total number of offsprings (i.e., 
impact degree) r is distributed according to the Borel distribution [21]: 

P{r) = (H r_1 V- ( 2 ) 



Using Stirling's formula (i.e., r! « y2jvrr r e r ), the above equation leads to 
the approximation 

P(r) oc r -3/V r(ln ^ +1) . (3) 



In particular, when \i — 1 (i.e., the critical case), the impact degree follows a 
power-law distribution with exponent —3/2. 



3.2.2 Power-law model 

In addition to the Poisson case, we consider the case where the number of 
offsprings is determined based on a power-law distribution (hereafter called 
Power-law model). Indeed, the number of offsprings for each reaction node 
is based on the outdegree in metabolic networks. Since real-world complex 



7 



networks including metabolic networks have power-law degree distributions 
[25|29f30j . the number of offsprings for each reaction node may also obey a 
power-law distribution. 

Saichev et al. [20|33] showed analytical asymptotic approximations for the 
distribution P(r) of the total number of offsprings (i.e., impact degree) r in 
the case where the number of offsprings d for each progenitor follows asymp- 
totically a power-law distribution P(d) oc l/d 1+1 and the mean number of 
offsprings has a given value \x. In particular they derive the following approx- 
imation for large impact degree r, in the case 1 < 7 < 2 where the variance 
of the number of offsprings is infinite: 

P( r ) / (l-/i)r-/i-l \ 

where 

oc 

<y9 7 (x) = J exp 


and v = /i(7 — l)7 1 / 7 " 1 r(— 7) 1 / 7 , where T(x) is the Gamma function. 

When fi = 1 (i.e., the critical case), in particular, the distribution of the total 
number of offsprings obeys the power-law distribution: P(r) oc l/r 1+1 / 7 . 

In addition, P(r) in the case of 7 > 2 is approximately similar to that in the 
case of Poisson model [i.e., Eq. (J2])] [20] . 

3.2.3 Empirical model 

The above probability distributions may be unsuitable to approximate real- 
world offspring distributions. In addition, we also consider a branching pro- 
cess model using empirical offspring distributions (hereafter called Empirical 
model). This way, we can estimate the distribution P(r) of the total number 
of offsprings (i.e., impact degree distributions) without the approximation of 
offspring distributions, although the model is analytically intractable. More- 
over, we can evaluate whether the prediction accuracy of the proposed method 
is influenced by the approximation of offspring distributions or the fidelity of 
branching processes. 

Let F(s) be the probability generating function of the impact degree r (i.e., the 
total number of offsprings), the function F(s) satisfies the recursive relation 



u 1 cos 



7T7 



COS 



7 ' f ^ J_ 

u ' sin I — + ux 



2 



du 



(5) 
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F(s) = f(sF(s)), 



(6) 



where f(s) denotes the probability generating function of the number d of 
offsprings of each node. 

Using the Lagrange expansion and the relation of P(r) = (l/r\)d r F(s)/ds r \ s= o, 
the distribution P(r) (i.e., impact degree distribution) is derived from the 
above implicit equation as the following explicit equation [20f2T] : 

(r > 0), (7) 

s=0 

where f(s) = X^=o x P(d)s d . The value <i max indicates the maximum of d, and 
the function P(d) corresponds to the probability density function of empirical 
d (i.e., empirical offspring distribution). In addition, P{r) = /(0) when r = 0. 

3.2.4 Parameter extraction 

To apply the Poisson and the power-law models, we need to estimate the 
model parameters, namely the mean \i and the exponent 7, from real metabolic 
networks. 

We estimate the mean of the number of offsprings for each reaction node by 
the empirical average: 

1 N 

^=JTrY, d ^ ( 8 ) 

where N is the total number of reaction nodes in a metabolic network. 

We estimate the exponent of a power-law offspring distribution using the max- 
imum likelihood estimation method [M] : 

-1 

(9) 



where N* is the set of reaction nodes with di > 0, and \N* \ indicates the total 
number of such reaction nodes. d m m is the minimum of di in the set of N*. 



P(r) 



1 d r 



r-l 



r! <is r_1 



f r (s 



df(s) 
ds 



7 



\N*\ 
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4 Evaluation of the branching process approximation 

We evaluated the above estimation methods for the impact degree distribu- 
tions on several real metabolic networks. 

We selected two bacteria [Escherichia coli (eco) and Bacillus subtilis (bsu)] 
and two eukaryotes [Saccharomyces cerevisiae (see) and Homo sapiens (hsa)] 
whose metabolic pathways have been well-identified. We downloaded the data 
of their metabolic networks, represented as bipartite networks as shown in Fig. 
[TJ from the KEGG database (version 0.7.1) [2|STj . The parenthetic three- letter 
codes correspond to KEGG organism identifiers [32] . 

Based on the KEGG metabolic network data, the impact degree distribu- 
tions in the metabolic networks were calculated using the method explained 
in Sec. |5J Moreover, we constructed the reaction networks of these species, 
and obtained the offspring distributions. Using Eqs. OH]) and OH]), the model 
parameters \i and 7 were extracted (see Table [3D- All metabolic networks show 
1 < 7 < 2, implying that the assumption of power-law model is suitable. 

Table 3 

Model parameters extracted from real metabolic networks. The character # indi- 
cates "the number of . The number of reaction nodes corresponds to the number 
of parents in branching processes. 



Species 


# Reaction nodes 


Mean u 


Exponent 7 


Escherichia coli 


1,467 


0.68 


1.52 


Bacillus subtilis 


1,279 


0.60 


1.65 


Saccharomyces cerevisiae 


1,172 


0.54 


1.73 


Homo sapiens 


1,982 


0.58 


1.50 



To test the validity of the offspring distribution models, we compared the 
fitting results on empirical offspring distributions in real metabolic networks 
between the power-law distribution and the Poisson distribution. Fig. [3] shows 
the cumulative offspring distributions from the real metabolic networks and 
the cumulative representation of theoretical distributions. The figure clearly 
indicates that the power-law distributions are more appropriate for modeling 
offspring distributions than the Poisson distributions. However, the power-law 
distributions may be not the best model because of the poor fittings for the 
larger d in the case of H. sapiens (Fig. EP) 

Using Eqs. ©, ©, and (J7J), we obtained the estimated impact degree distribu- 
tions using the branching process approximation. Fig. [4] shows the comparison 
between the observed cumulative impact degree distributions and estimated 
ones. The theoretical predictions (lines) are in good agreement with the real 
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Fig. 3. Cumulative offspring distributions of Escherichia coli (A), Bacillus subtilis 
(B), Saccharomyces cerevisiae (C), and Homo sapiens (D). P cum (x) is defined as 
P(X > x). Note that P cum (l) < 1 because the cumulative distributions also con- 
sider the case of d = 0. P CU m(_0) is not shown due to the logarithmic display. The 
symbols indicate observed data. The black solid lines and dashed lines correspond 
to the cumulative representations of the power-law distribution with the exponent 
estimated by Eq. ([9]) and the Poisson distribution with the mean obtained from Eq. 
([8]), respectively. 

impact degree distributions (symbols), suggesting the relevance of branching 
process approximations. Note that the impact degree distributions does not 
follow a clear power law and show an exponential cut-off for larger impact 
degrees because ji < 1 (i.e., not the critical case). Eqs. (j3J) or PJ can explain 
this distributional tendency. 

To evaluate the prediction accuracy for impact degree distributions between 
the Power-law model and the Poisson model, we evaluated the distributional 
distance (i.e., Kolmogorov-Smirnov statistics) between the observed distribu- 
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Fig. 4. Cumulative impact degree r distributions of Escherichia coli (A), Bacillus 
subtilis (B), Saccharomyces cerevisiae (C), and Homo sapiens (D). The symbols 
indicate observed data. The black solid lines and dashed lines correspond to the 
cumulative representations of theoretical distributions of the power-law model and 
the Poisson model, respectively. The gray solid lines are the the cumulative repre- 
sentations of theoretical distributions from the empirical model. Note that P cum {x) 
is defined as P(X > x). 



tions and theoretical distributions (Table HJ). The power-law model is better 
than the Poisson model on all networks. The empirical model outperforms 
the power law model on both bacterial networks. Surprisingly, the power law 
model outperforms the empirical model on both eukaryote's networks. 
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Table 4 

Prediction accuracy for the impact degree distributions: Kolmogorov-Smirnov (KS) 
distance, defined as sup,,, \ R(x) — M(x)\, where R(x) and M(x) are empirical distri- 
butions and theoretical distributions, respectively. The parenthetic values indicate 
the logarithmic P- values p from the KS test, defined as — log 10 (p). The emphasized 
values correspond to the best accuracy. 



Species 


Poisson model 


Power-law model 


Empirical model 


Escherichia coli 


0.11 (2.06) 


0.08 (0.86) 


0.05 (0.27) 


Bacillus subtilis 


0.18 (4.65) 


0.08 (0.75) 


0.05 (0.15) 


Saccharomyces cerevisiae 


0.14 (2.76) 


0.07 (0.40) 


0.07 (0.54) 


Homo sapiens 


0.12 (3.12) 


0.09 (1.86) 


0.14 (5.09) 



5 Discussion and conclusion 



We proposed a model to estimate the impact degree distributions in metabolic 
networks, using a branching process approximation, and demonstrated its va- 
lidity on real data. 

The power-law model could more accurately estimate the impact degree distri- 
butions in real metabolic networks than the Poisson model because the number 
of offsprings for each reaction node is assumed to follow the power-law dis- 
tribution. Especially, the power-law model showed the significant agreements 
between the predicted distribution and the observed distribution although the 
case of H. sapiens represented the small P- value for the KS test (i.e., the low 
probability that the distribution is similar between models and observed data). 

However, there is no great difference of the prediction accuracy for estimating 
the impact degree distributions between the power-law model and the Poisson 
model; thus, the Poisson model may be useful for a rough estimate of the 
impact degree distributions. 

Intrinsically, the distribution of the total number of offsprings [i.e., P(r)] is not 
significantly different between the power-law model and the Poisson model in 
the case of smaller r. As a simple example, we here consider the case of \x = 1 
(i.e., the critical case). In this case, the impact degree distributions P(r) of 
the power-law model and the Poisson model correspond to oc r~ < - 1+1 / 7 ^ and 
oc r -3 / 2 , respectively. Especially, P{r) of the power-law model is the power- 
law distribution with the exponent ranging between —2 and —1.5 because of 
1 < 7 < 2; thus it is not critically different from r~ 3 / 2 of the Possion model for 
smaller r. Since the impact degrees of real metabolic networks (i.e., r) were 
relatively small (r < 60), there might be no great difference of the prediction 
accuracy between these models. 
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The similarity of predicted distributions of the total number of offsprings be- 
tween the Poisson model and the power-law model is also explained using the 
Otter's theorem [Theorem 13.1 in [21] ]. This theorem indicates that the dis- 
tribution P(r) of the total number of offsprings has the universal property of 
a power-law tail with the exponent —3/2 as r -)■ oo under mild conditions on 
offspring distributions, and it implies that the types of offspring distributions 
hardly influence the distribution of the total number of offsprings (i.e., impact 
degree distributions). Note the Otter's theorem does not contradict with the 
analytical distribution P(r) of the power-law model because this theorem is 
not directly applicable to the power-law model due to the different assump- 
tions in the derivation of P(r) between the power-law model and the Otter's 
theorem. 

The prediction performance is influenced by the assumption of offspring distri- 
butions and the fidelity of branching processes. To purely evaluate the validity 
of branching process approximation, the empirical model is useful because of 
no assumption of offspring distributions. It is expected that the empirical 
model show the best prediction accuracy because of using empirical offspring 
distributions. In the case of bacteria (i.e., E. coli and B. subtilis), this expec- 
tation is true, suggesting that the branching process approximation is useful 
for estimating the impact degree distributions. In the case of eukaryotes (i.e., 
S. cerevisiae and H. sapiens), on the other hand, we observed the unexpected 
results: the empirical model shows the relatively-low prediction accuracy. Es- 
pecially, the prediction accuracy of the empirical model is lowest in the case 
of H. sapiens (human). This result implies limitations to the estimation of 
impact degree distributions based on the branching process approximation in 
the case of metabolic networks of eukaryotes (i.e., higher organisms). 

A limitation of the branching process approximation is that we need to assume 
tree structures of networks (i.e., no cycles). The presence of cycles may lead 
to an overestimation of the number of offsprings di for each reaction, because 
some offsprings of a progenitor (reaction node) may have already been inacti- 
vated due to cycle structures. From this reason, the models may overestimate 
the impact degree distributions. On the other hand, however, some reactions 
with more than one incoming edge may be inactivated in the presence of cy- 
cles, if all their parents are inactivated through different paths. In this case, 
the number of offsprings is underestimated. The estimation of the number 
of offsprings depends on the relative importance of these two effects, and it 
may be not simple. The difficulty in this estimation is also a limitation of the 
model. 

The branching process model needs to be improved by considering additional 
assumptions other than metabolic network structures in order to obtain bet- 
ter predictions. For example, the assumption of variable propagations in the 
branching process [33], in which the mean of offspring distributions differs 



14 



at each propagation stage, may be useful because the degree of propagation 
may depend on metabolic dynamics such as gene expressions and metabolite 
concentrations. To apply this modified branching process model, time-series 
data on metabolic dynamics after gene disruptions, which are obtained by 
metabolomic analysis, are necessary for estimating the mean of offspring distri- 
bution at each propagation stage (i.e., time after gene disruptions). Since such 
data are unavailable at present, however, it is difficult to apply this modified 
branching process model. Similarly, other biologically-suitable assumptions are 
hardly determined because of few observed data on metabolic dynamics. Al- 
though there are above constraints on observed data on metabolic dynamics, 
we believe that the consideration of the additional information improves the 
prediction of impact degree distributions. In the future, the improvement of 
the prediction of the impact degree distributions using the branching process 
approximation may be possible with available data on metabolic dynamics. 

The branching process approximation is useful for estimating the impact de- 
gree distributions in metabolic networks although it has the above limitations; 
thus, it may be a powerful tool for evaluating a robustness of biological sys- 
tems. 
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